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A program  is  given  that  numerically  inverts  a given  Laplace 
transform  for  discrete  times  specified  by  the  user,  and  (most 
of  the  time)  to  within  a specified  absolute  error.  The 
computation  is  adaptive  in  that  the  program  decides  whether  to 
compute  a given  point  by  integrating  the  Bromwich  integral 
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f(s)  ds,  t^O 


or  to  interpolate  the  point  from  the  set  of  those  already 
computed  by  integrating  the  Bromwich  integral.  If  the 
integration  is  performed,  it  is  done  by  accelerating  partial 
sums  to  the  limit,  with  the  partial  suras  obtained  by  Gaussian 
quadrature  with  error  control.  Three  examples,  including  two 
in  which  the  transform  is  not  the  quotient  of  polynomials, 
indicate  that  the  reliability  of  the  program  is  good. 
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1.  INTRODUCTION 

The  Laplace  transform  is  used  frequently  in  engineering  amalysis 
primarily  because  many  problems  are  easily  formulated  with  it.  The 
final  solution  requires  inverting  a tranform  into  ,the  time  domain, 
sometimes  not  possible  in  terms  of  common  functions.  Hence,  many 
numerical  inversion  techniques  have  evolved.  Most  of  the 
techniques  use  the  basic  definition 


= r 


e~®^  F(t)  dt,  t>0  , 


evaluate  f(s)  for  selected  real  s.  and  specific  tj,  and  set  up  a matrix 
equation  for  solution,  perhaps  after  a transformation  to  allow  the  use 
of  orthogonal  functions  to  simplify  the  solution.  While  these  methods 
are  extremely  fast  when  done  on  a computer,  they  are  limited  in  the 
accuracy  obtainable  and  the  class  of  f (s)  to  which  they  can  be  applied. 
Functions  of  time  that  are  discontinuous  or  "nonsmooth"  in  some  manner 
cannot  be  readily  obtained.  Berger  and  Duangudom  also  point  out^  that 
Berger's  method^  may  even  suffer  from  accuracy  loss  both  for  large  t and 
smooth  oscillatory  functions  such  as  F(t)  = sin  t. 

The  techniques  referenced  above  are  not  applicable  to  Laplace 
transforms  obtained  in  the  emalysls  of  certain  fluid  transmission  lines, 
where  they  are  of  the  form 


f(s)  = 


s cosh  lG(s)l 


The  inverse  of  equation  (2)  may  be  discontinuous  [depending  on  the 
nature  of  G(s)]  because  of  wave  reflections  in  the  line.  Although  these 
problems  can  also  be  formulated  for  solution  by  the  method  of 


^C.  Lanczos,  Applied  Analysis,  Prentice-Hall,  Inc.,  Englewood  Cliffs, 
NJ  (1956) pp.  284-303. 

^A.  Papoulis,  "A  New  Method  of  Inversion  of  the  Laplace  Transform," 
Quarterly  of  Applied  Mathematics,  1£  (1957). 

^3.  S.  Berger,  "Inversion  of  the  N-Dimensional  Laplace  Transform," 
Mathematics  of  Computation,  20  (1966),  pp.  418-421. 

‘*R.  Bellman,  R,  E.  Kalaba,  and  J.  A.  Loclcett,  Numerical  Inversion  of 
the  Laplace  Transform,  American  Elsevier  Publishing  Co. , New  York 
(1966) . 

^B.  S.  Berger  and  S.  Duangudom,  “A  Technique  for  Increasing  the 
Accuracy  of  the  Numerical  Inversion  of  the  Laplace  Transform  with 
Applications ,"  ASME  Journal  of  Applied  Mechanics,  paper  No.  73-WA/APM-l 
(December  1973). 
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characteristics  and  other  methods,  it  was  felt  that  a general-purpose 
inversion  progreun  for  arbitrary  Laplace  transforms  would  have  excellent 
utility. 

The  method  developed  is  similar  to  that  of  Schmittroth^  except  that 
an  error  control  is  provided.  It  is  based  on  integration  of  the 
Bromwich  integral  formula^ 


F(t) 


2iri  J 

Y-i“> 


st 


e 


f(s)  ds,  t?<0  , 


(3) 


where  s = x + iy  (x,y  real)  and  the  integration  in  the  complex  plane  is 
performed  along  the  line  x = y«  The  number  Y is  arbitrary  but  must  be 
chosen  so  that  the  line  x = y lies  to  the  right  of  all  singularities 
(poles,  branch  points,  or  essential  singularities). 

This  method  is  slower  than  other  methods  because  F(t)  is  obtained 
point  by  point  by  computing  eui  infinite  integral.  To  help  speed  matters 
along,  the  final  algorithm  contains  an  adaptive  scheme  that  decides 
vrtiether  or  not  to  interpolate  a given  point  from  the  set  of  those  then 
existing.  Despite  the  relative  slowness,  good  results  have  been 
obtained  because  of  the  current  high  speeds  of  digital  computers.  Ease 
of  use  and  general  utility  of  the  program  are  its  prime  assets. 


2.  TRANSFORMATIONS  TO  A REAL  INTEGRAL 
With 


s = Y + iy  ^4) 

ds  = i dy  (5) 


equation  (3)  transforms  to 


I 


M 


F(t) 


2ir 


iyt 


f(Y+iy)  dy  , 


(6) 


®L.  A.  Schmittroth,  "Numerical  Inversion  of  Laplace  Transforms,"  ACM 
Communications  (March  1960),  pp.  171-179. 

^F.  Scheid,  Theory  and  Problems  of  Numerical  Analysis,  McGraw-Hill 
Book  Co.,  Shaum's  Outline  Series,  New  York  (1968),  p.  125  ff. 


(cos  yt  + i sin  yt)  f(>+iy)  dy 


With  the  transformation  y = -y'  in  the  first  integral,  (7)  reduces  to 


{ (cos  yt  + i sin  yt)  f (y+iy) 


We  now  make  use  of  conjugate  properties  of  complex  numbers 


If  f(s)  is  expressed  by  a Taylor  series  (we  assume  analyticity) 


we  have 


t*! ' 


Further  assuming  the  a. 


real ,® 


then 


[f(s)l*  « a.(s*)^  = f(s*) 

j=0  ^ 


(12) 


writing 


f(y+iy)  = Re(f)  + i Im(f) 


(13) 


we  have  from  equation  (12) 


f (y-iy)  = Re(f)  - i Im(f) 


(14) 


Equation  (8)  reduces  to 


e^^  r 

F(t)  = — — / {Re(f(Y+iy)l  cos  yt  - Im[f(Y+ly)]  sin  yt} 

7i 


dy 


(15) 


Since  F(-t)  = 0,  the  component  parts  of  equation  (15)  are  equal,  but  of 
opposite  sign.  Thus 


2e^^  f 

F(t)  = I Re[f(Y+iy)]  cos  yt  dy 

TT  J 


(16) 


2^/ 

IT  -o. 


lm[f(Y+iy)]  sin  yt  dy 


We  use  the  first  integral  in  equation  (16)  as  the  basic  one  to  be 
evaluated.  With  the  final  transformation 


0)  = yt 


do)  = t dy  , 


(17) 


®R.  V.  Churchill,  Complex  Variables  and  Applications,  2nd  ed., 
McGraw-Hill  Book  Co.,  New  York  (1960).  (The  result,  eg  (12),  is  the 
so-called  "Principle  of  Reflection."  The  simplest  test  to  determine  if 
a is  real  is  that  f(s)  is  real  whenever  s is  real.) 


i 


i 


i 

'4 

1 

I 


I 


our  transformed  integral  becomes 


F(t)  = 


f-K-a] 


cos  b)  du 


3.  COMPUTING  THE  TRANSFORMED  INTEGRAL 

Subroutine  TPOINT  computes  the  transformed  integral  [eq  (18)]  (see 
the  listing  in  appendix  A for  a description  of  the  argiunents) . To  do 
this,  the  infinite  integral  is  first  changed  to  an  infinite  sum  of 
finite  integrals.  Arbitrarily  taking  one  cycle  of  the  cos  w factor  as 
the  range  of  integration  for  each  finite  integral,  we  define 


so  that 


1 

^ ‘211  (j-1) 


a] 


cos  0)  do) 


F(t)  = 


E F (t)  . 

j=l  ^ 


Two  problems  remain  in  computing  F(t).  First,  each  Fj(t)  must  be 
computed  accurately.  Second,  the  infinite  sum  in  equation  (20)  must  be 
changed  to  a finite  sum,  say  over  N terms.  To  minimize  N,  we  also 
intend  to  apply  some  nonlinear  transformation  algorithm  to  the  sequence 


X of  partial  sums, 
m 


m 

= Z 

j=i  ^ 


in  order  to  accelerate  it  to  the  limit  as  m-+«>. 


3.1  Computing  F^  (t) 

After  some  experimentation  with  QATR,  thw  IBM  360  scientific 
routine^  to  perform  quadrature  integration,  it  was  decided  that  the 
routine  was  inadequate  to  the  task  of  computing  (t) . To  obtain  fast 

^"Systmi/360  Scientific  Subroutine  Package,  Version  III,  Programmer’s 
Manual,"  Application  Program  GH20-0205-4  (1968),  pp.  297-298. 
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convergence  in  computing  equation  (20) , it  is  desirable  to  compute 
with  Y close  to  the  singularities.  This,  in  turn,  produces  large  peaks 
in  the  integrand  of  equation  (19)  so  that  QATR  frequently  returned  with 
the  error  code  indicating  that  the  accuracy  could  not  be  reached  because 
of  rounding  errors.  In  that  case  the  smallest  distance  of  y from  the 
singularities  was  doubled  and  the  procedure  repeated.  The  procedure 
turned  out  too  time-consuming  even  for  the  simplest  transforms. 


The  following  procedure  was  adopted  to  compute  equation  (19) 
for  each  j ; 


(1)  Compute  equation  (19)  with  an  m^-point  Gaussian  quadrature 
formula,  vAiere  m^^  is  either  6,  12,  16,  24,  32,  40,  48,  64,  80,  or  96, 
for  k = 1,  2,  ...,  10,  arbitrarily  using  = 32  for  j = 1.  Compute 
again,  using  an  formula.  If  the  two  results  are  and  we 
require  that 


(?-.) 


where  E is  an  absolute  error  parameter  supplied  by  the  user. 


(2)  If  equation  (22)  is  satisfied,  (t)  is  taken  to  be 
If  not,  the  svibscript  of  m is  increased  until  equation  (22)  is 


satisfied. 


(3)  The  first  pair  of  subscripts  for  which  equation  (22)  is 
satisfied  is  also  used  in  the  next  interval  (j  increased) . If  in 
addition 


S - Vi  i ■ 


(23) 


the  subscript  k is  decreased  by  one  for  the  initial  try  in  the  next 
interval  (unless  k = 1) . 

(4)  If  k = 9 and  equation  (22)  is  not  satisfied,  the  distance 
that  Y is  from  a singularity  is  dov±)led  and  step  1 is  repeated. 

(5)  If  Y t >11,  the  procedure  is  aborted.  (This  has  never 
happened . ) 


TO  confute  equation  (19)  by  a Gaussian  quadrature  formula,  we  j 

first  let 

(24) 

Thus, 


(25) 

»=  —j  Re  ^f  3-  ^^^jcos  W2  dz  , 


since  2j-l  is  an  odd 

7 

computes 


I where  the  abscissas  z and  the  weights  W*  are  well  tabulated  for  a 

^ n n 

[ variety  of  m^^.  It  is  convenient  to  store  the  data 

( 

( P - irz  (27) 

n n 


integer.  Xn  Gaussian  quadrature  formula 

I - g(z)  dz  s.  £ g(z„)  , (26) 

n-1 


u - S(z+1)  ♦ 2s(j-l) 


to  transform  the  limits  of  z from  12s(j-l),  2irjl  to  1-1,1). 
equation  (19)  tramsforms  to 

F^(t)  « J Re  ^f  cos  (ir(z+2j-l)]  dz 


to  eliminate  this  multiplication  during  the  course  of  the  solution. 
Similarly,  it  is  convenient  to  store 

W » - {cos  P ^ W*  (28) 

n V n/  n 

since  these  are  also  independent  of  j.  From  equations  (25)  to  (28),  the 
Gaussian  quadrature  formula  used  is  then 


V.  Scheid,  Theory  and  Problems  of  Numerical  Analysis,  McGraw-Hill 
Book  Co.,  Shaum’s  Outline  Series,  New  York  (1968)  , p.  125  ff. 
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(29) 


>"k  r / ijp  +TT{2j-i)i\i 

yt,  . E w„  Ke  [f  (r* 


The  algorithm  described  is  very  efficient,  once  the  correct  Ic 
is  found  in  the  first  interval,  because  success i/e  intervals  are  usually 
similar  to  one  another,  and  there  is  never  any  need  to  compute  the 
cosine  function. 

Condition  (22)  does  not  imply  that  is  computed  to  vithin  an 
error  E,  but  experience  has  shown  that  is  the  case  for  most  of  the 
exaii^>les  tried.  Errors,  however,  can  accumulate  in  X^.  The  final  error 
in  F(t)  is  subject  to  these  errors  and  the  errors  due  to  the 
acceleration  methods  described  below. 

3.2  Accelerating  the  Sequence  X to  the  Limit 

n 

The  sequence  X in  equation  (21)  obtained  from  summing  F.(t), 
m 3 

j = 1,  2,...,  m,  approaches  the  limit  F(t)  as  m-xo,  if  y is  to  the  right 

of  all  singularities.  This  is  the  result  of  the  Bromwich  integral 

theorem.  Since  convergence  may  be  very  slow  it  is  prudent  to  consider 

some  acceleration  technique.  Svibroutine  TEAS^  seemed  ideal  for  this 

purpose.  This  subroutine,  however,  was  frequently  fooled  by  the 

sequences  obtained  from  transforms  of  the  type  given  in  equation  (2) . 

After  some  experimentation  with  the  Shanks  algorithm,^®  upon  which  TEAS 

is  based,  it  was  decided  that  the  algorithm  was  too  sensitive  to 

triincation  errors,  since  each  X -X  , could  have  an  absolute  error  of 

m m-1 

cdx>ut  E.  It  was  therefore  decided  to  simply  apply  an 

ej  transformation (i.e.,  an  Aitken  6^  transformation)  to  a modified 

svibsequence  of  X , as  descrilaed  below, 
n 


’^"System/ 360  Scientific  Subroutine  Package,  Version  III,  Programmer's 
Manual,"  Application  Program  GH20-0205-4  (1968),  pp.  234-237. 

Shanks,  "Son-Linear  Transformation  of  Divergent  and  Slowly 
Convergent  Sequences,"  2£  (1955),  pp.  1-42. 
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m = 


Aitken's  6^  transformation  generates  a sequence 

3,  4,...,  N,  from  a sequence  X , m = 1,  2,  with 

m 


m 


2m  = 


m 


fX  -X 
V m m-1/ 

X -2X  ,+X  - 

m m-1  m-2 


(30) 


If  m«>  X “ L,  and  L-x  approaches  zero  nearly  geometrically,  then  0 
in  ro  iR 

approaches  the  limit  L faster^®  than  X . This  nonlinear  transformation 

m 

is  frequently  used  to  accelerate  infinite  sums  and  can  be  extremely 
effective  to  this  end.  The  subsequence  to  be  chosen  is  aimed  at 
obtaining  differences  that  are  nearly  geometric. 

Figure  1 shows  three  types  of  sequences  obtained  by  the 
procedvure  described.  Type  (a)  eventually  converges  monotonically  to  a 
limit  with  the  absolute  difference  between  successive  values  becoming 
smaller.  Type  (b)  is  similar  to  type  (a)  except  that  the  limit  is 
approached  via  oscillation.  Type  (c)  is  a combination  of  (a)  and  (b) , 
but  the  oscillation  is  such  that  the  maximums  get  larger  or  the  minimums 
get  smaller. 

More  complicated  oscillatory  sequences  are  obtained  with  some 
Laplace  transforms,  and  it  is  difficult  to  devise  an  algorithm  that  will 
never  be  fooled.  For  excunple,  figure  2 shows  a sample  sequence  actually 
obtained  for  equation  (2)  where  G(s)  = /s+s^ . After  a small  niimber  of 
intervals,  it  appears  to  be  a type  (b)  sequence,  but  actually  the 
high-frequency  oscillation  is  superimposed  on  a low-frequency 
oscillation,  and  the  limit  is  considerably  greater  than  the  apparent 
one. 


In  order  to  minimize  the  false  limit  possibility,  the  following 
procedure  is  used; 

(1)  A minimum  of  11  points  is  used,  i.e.,  the  integration  is 
carried  to  at  least  22 tt. 

(2)  A type  (a)  sequence  is  assumed  if  11  consecutive  points  are 
roonotonic  with  successive  absolute  differences  becoming  smaller. 


Shanks,  "Non-Linear 
Convergent  Sequences 34 


Trans formaticn  of 
(1955),  pp.  1-42. 


Divergent 


and 


Slowly 
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limit 


(a)  Eventually  approaches  a Irnut  monotonically  (successive 
diffezences  get  smaller) . 


(b)  Eventually  approaches  a limit  via  oscillation  where  the 
envelopes  of  the  maximums  get  smaller  and  the  minimums 
get  larger. 


• - • • 


(c)  Eventually  approaches  a limit  via  oscillation  where  the 
envelopes  of  the  maximums  get  larger  or  the  minimums  get 
smaller. 

Figure  1.  Typical  sequences  X^. 


Figure  2.  Example  of  a complicated  sequence. 


(3)  Peaks  (extremums)  are  stored  in  a peak  array  P.  A peak  is 
assumed  if  2^)  (^i-l”^i-2)  ^ pareQsola  is  passed  through  the 
three  points  X^,  and  X^_2  to  determine  a corrected  peak.  If  P(J) 
is  a new  peak  and  it  is  determined  that  it  is  a maximum  (or  minimum) , 
its  value  must  be  less  (or  greater)  than  P(J-2)  to  be  placed  in  the  P 
array . 


(4)  otherwise  P(J)  replaces  P(J-2)  and  P(J-2)  is  stored  in  an 
envelope  array  R(K) . A type  (b)  sequence  for  the  P array  is  assumed  if 
J reaches  at  least  5.  A type  (c)  sequence  is  assumed  if  K reaches  11 
and  consecutive  R values  are  monotonic  with  successive  absolute 
differences  beccxning  smaller. 

(5)  No  matter  what  type  of  sequence  is  finally  assumed,  the 
<S^  transformation  is  applied  to  the  last  five  values  of  the  X,  P, 
or  R sequence.  The  resulting  three  extrapolation  values  must  all  be 
within  E of  each  other  for  the  result  to  be  accepted.  Otherwise  more  X. 
are  computed.  ^ 

(6)  If  the  accepted  extrapolated  values  oscillate,  the 
6^  transformation  is  applied  to  them  to  get  a final  answer.  Otherwise, 
the  last  extrapolated  value  is  taken  as  the  answer. 

Although  the  above  procedure  is  not  foolproof,  it  has  rarely 
failed  to  give  a satisfactory  answer,  because  the  niagnitude  of 
differences  of  the  final  array  used  in  the  transformation  gets 
smaller,  i.e.,  is  nearly  geometric. 

4.  COMPUTING  MANY  F(t)  WITH  SUBROUTINE  POINTS 

Subroutine  POINTS  computes  F(t)  at  all  t supplied  in  the  array  Tl 
(see  listing).  The  Tl(l)  should  be  ordered  so  that  the  values  increase 
with  I (for  I < 501) . 

To  save  computer  time,  provision  is  made  to  compute  some  of  the 
answers  Y1(I)  by  interpolation.  The  user  supplies  the  number  N1  of 
total  points  (the  dimension  of  Tl)  and  the  minimum  number  N2(>2)  that  he 
desires  to  be  computed  by  integrating  the  Bromwich  integral  as  described 
in  section  2.  POINTS  distributes  these  as  equally  spaced  as  possible. 
From  these  computed  Y1 (I) , it  then  attempts  interpolation  for  all  other 
points  using  subroutine  INTERP,  based  on  the  routine  ALI  in  the  360 


scientific  subroutine  package.^  The  modifications  producing  INTERP 
allow  for  the  natxiral  ordering  of  the  T1  array  and  do  not  destroy  that 
array.  The  listing  is  fully  documented;  changes  of  ALZ  are  identified 
by  the  lack  of  identification  in  columns  73  to  80. 

If  INTERP  cannot  obtain  an  accurate  enough  interpolated  value  (using 
the  absolute  error  requirement  specified  by  the  user)  for  any  reason, 
the  point  is  calculated  by  calling  TPOINT.  Thus,  a dense  array  T1  can 
be  specified  (say  for  plotting  purposes) , but  only  as  many  points  as 
will  be  required  will  be  calculated  by  integration. 

Provision  is  also  made  to  enter  data  (as  is  required  when  t < 0)  and 
to  force  integration  of  the  Bromwich  integral  for  any  point,  by  means  of 
the  I/O  array  IE.  See  the  listing  for  a description  of  its  use.  Thus, 
suspect  points  computed  by  interpolation  can  be  recomputed. 

Another  provision  is  made  for  printing  the  results  either  as  they 
are  computed  (not  in  the  order  of  increasing  Tl)  or  in  ordered  form 
after  all  confutations  are  finished.  Printing  can  also  be  suppressed. 
The  varicd^le  IPRINT  controls  these  features. 


5.  EXAMPLES 

Three  examples  are  given  to  demonstrate  the  use  of  POINTS. 
5.1  Example  1 

We  do  first  a trivial  example 


f(s) 


s**  - 6s^  + 1 

(s2  + !)*♦ 


(31) 


the  solution  of  which  is  F(t)  = t^  (cos  t)/6.  We  wish  the  func  .onal 
values  for  the  range  0 < t < 10.  Since  the  maximum  magnitude  of  the 
function  is  edx>ut  200  in  this  remge,  a good  graph  of  the  function  is 
obtained  with  an  absolute  error  specification  of  0.05.  About  200  points 
provide  a dense  enough  grid  for  this  case  (we  use  201  to  get  At  = 0.05) 
and  we  limit  the  number  of  intervals  to  integrate  to  100  per  point 
because  F(t)  is  so  smooth.  A main  program  to  print  points  as  they  are 
calculated  (IPRINT  = 1)  is  shown. 


^"System/360  Scientific  Subroutine  Package,  Version  III,  Programmer's 
Manual,"  Application  Program  GH20-0205-4  (1968),  pp.  241-242, 
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DIMENSION  Wl(lOO),  T(201) , Y(201),  IE(201) 


5 T(I)  - FLOAT  (I-D/20 


(a)  The  IE  and  T eirrays  must  be  defined  before  POINTS  Is  called; 
IE  (I)  s 0 in^lles  that  T(I<  must  l>e  calculated. 


(b)  Since  T(l)  = 0,  Y(l)  must  be  supplied.  This  can  usually  be 
found  easily  from  the  initial  value  theorem  [F(0>  = sf(s)].  Y(l)  is 

set  to  zero,  amd  IE(1)  s i to  indicate  that  the  data  are  supplied. 


(c)  N2  = 21  was  chosen  arbitrarily  as  the  number  of  points  to  be 
computed  by  integration  initially.  (Of  course,  if  supplied  with 
IE  (I)  »=  1 the  integration  is  bypassed.) 


Running  times  will  be  proportional  to  the  time  of  computing 
f (s) , so  any  procedure  to  speed  the  process  is  helpful . The  Laplace 
transform  routine  T1  was  written  with  this  in  mind: 


SUBROUTINE  Tl(S,F) 
COMPLEX  S,  F,  S2,  SI 
S2  = S*S 


SI  = S2  + 1. 

SI  = S1*S1 

SI  = S1*S1 

F = (1.  + S2*(S2  - 6.))/Sl 

RETURN 

END 

Figure  3 (p.l9)  shows  a partial  output.  Point  numbers  1 

through  201  in  steps  of  10  (uncaptioned  column  at  left)  are  computed 
first.  Point  1 was  supplied  and  not  computed  by  integration  so  that 
IE(1)  was  returned  equal  to  1.  The  other  points  were  computed  by 
integration  since  IE(1)  > 11.  The  fact  that  IE(1)  was  11  or  12  shows 
that  convergence  was  very  fast. 

POINTS  then  tries  to  interpolate  other  points  and  succeeds 
except  for  points  146  and  156,  which  are  computed  next  by  integration. 
The  program  then  reinterpolates  (using  the  new  data  also)  and  starts 
listing  the  remaining  points  with  IE  = 5.  The  entire  process  took  about 
2 s on  an  ia<  360/195  computer. 

5.2  Example  2 

We  invert  equation  (2)  with  G(s)  =Vs^  + s.  The  main  program 
is  similar  to  example  1,  except  that  the  call  to  POINTS  will  be 

CT^L  P01NTS(TRANS,0., .001,300,Wl,201,41,T,y,lE,2). 


Because  of  knowledge  of  the  system,  discontinuities  are  expected  and  41 
initial  points  are  requested  by  integration,  each  point  being  allowed  to 
sum  and  extrapolate  300  intervals.  The  error  parameter  w^s  chosen  to  be 
0.001  to  obtain  graphical  accviracy,  since  it  is  known  that  the  steady 
state  value  is  one.  IPRINT  was  set  to  2 to  list  all  values  in  sequence. 
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V4^ 


riME 


FUNCTION  EFftCT  CCCE 


1 0.0 

1 1 5.0C00C0CE-0  1 
21  l.OCCCCCCE  CC 
31  l.itCCCCCCE  00 
41  2.0CCCCCCE  CC 
5i  2.5CCCCCCE  CC 
61  3.CCCCCCCE  CC 
71  3.SCOOCOOE  OC 
81  >t.CCCCCCCE  CC 
91  A.bCCCCGCE  CC 

lOl  5.CCCCCCCE  CC 
UL  S.5C0CCCCE  CC 
121  6.CCC0CCCL  OC 
131  6.5CC0CCCE  CC 
141  7.0300C0CL  OC 

151  7.5CCCCCCE  CC 
161  8.0CCCCCCE  CC 
171  8.5CC0CCCE  CC 
181  9.CCCCCCCE  CC 
191  9.5CCCC0CE  OC 
2C1  l.CCCOCCCE  01 
146  7.25CCCCCE  CC 
156  7.75CCCCCE  CC 

2 4 .9999997E-G2 

3 9.9999964E-02 

4 1,4999S98E-C1 

5 1.9999999E-CI 

6 2.5CCCCCCE-CI 

7 2.9999995E-C1 

8 3.4S99996E-C1 

9 3.999999bt-01 

13  4.4S99999E-01 

12  5.4999995E-C1 

13  5.999999a-Cl 

14  t.499999EE-01 

15  t.9959999E-31 

16  7.5CCCCCCE-C1 

17  7.9999995E-01 


c.c 

1.82E2734E-Ci 

1 

11 

9.3373-37E-C2 

11 

i .98C5E59E-C2 

11 

-5.54E8718E-C1 

11 

-2.066375rE  CC 

11 

-4.4551315E  CC 

11 

-6.692C481E  CC 

11 

-6.9725543E  OC 

u; 

-2.2C58296E  CC 

11 

5.9134254E  CC 

12 

l.9t74'93E  Cl 

12 

3.436929CE  01 

11 

4.47C2C72E  Cl 

12 

4.3104736E  Cl 

12 

2.4276C07E  Cl 

12 

-1.2417i6t£  Cl 

12 

-6.1632935E  31 

12 

-1.1C72433E  C2 

12 

-1.425U85E  C2 

12 

-1.398759CE  C2 

12 

3.6C76E28E  Cl 

12 

e.0533C57E  CC 

12 

-5. 7944469E-C4 

5 

-6.22e3944£-C4 

5 

-1.3219^C3E-C4 

5 

S.925U5CL-C4 

5 

2.45  32527E-Ci 

5 

4,549C488E-C3 

5 

7.1798936E-C3 

5' 

1.3345794E-C2 

5 

1.4C46736E-02 

5 

2.3C53756E-C2 

5 

i.b359149£-C2 

5 

2.42CC989E-C2 

5 

4.357717CE-C2 

5 

4.748S41CE-02 

5 

7.6CC1465E-C2 

5 

Figure  3.  Partial  output  of  example  1. 
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The  transform  subroutine  is  trivial  to  vnrite: 
SOBRODTXNE  TRANS (S,P) 

COMPLEX  S,F,E 
E = CEXP(CSQRT(S*S  + S) ) 

F = 2./(S*(E  + l./E)) 

RETURN 

END 


Figure  4 shows  a partial  output  with  the  point  numbers  listed 
in  the  proper  sequence.  All  points  that  were  computed  by  integration 
were  done  successfully  (in  less  than  300  intervals)  with  CPU  time 
totalling  about  17  s on  cin  IBM  360/195.  The  maximvon  number  of  intervals 
required  was  284  for  point  103  (t  = 5.1) . Figure  5 shows  a graph  of  the 
x>mputed  function.  There  actually  are  discontinuities  at  every  odd 
integer  value  of  time,  6md  the  adaptive  procedure  clusters  mcuiy 
calculations  about  1,  3,  and  5.  The  interpolated  values  at  t = 0.95  euid 
t = 4.95  were  accepted  by  POINTS  and  are  in  error.  Similarly,  points 
near  t = 7 ^md  t = 9 were  computed  by  interpolation  eUid  are  in  error. 
Any  suspect  points  can  be  forced  to  be  computed  by  integration  with  the 
input  IE  * 3 (see  pp.  21  and  22  for  fig.  4 and  5) . 

5.3  Example  3 


As  a final  example,  we  Invert  a much  more  complex  but 
physically  realistic  case  for  equation  (2),  with 


G(s)  = s 


\l  1 \ 

/ .8Jj{vC73,F)  V 

_ 2Ji(F)  J 

FJo (F)  ! 

\ vCt1fJo(  /TTit)/ 

where 


F = i/is,  i = , 


(32) 


(33) 


and  Jo (A)  and  Ji (A)  are  Bessel  functions  of  the  complex  argvonent  A. 
G(s)  can  be  shown  to  be  real  when  s is  real;  it  is  a legitimate 
transform  with  a real  inverse. 
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TIME 


FUNCTION  ERROR  CODE 


1 

0.0 

0.0 

1 

2 

4.9S99997E-02 

2.3233570E-04 

5 

3 

9.9999964E-32 

4.9645780E-04 

5 

4 

1 .4999998E-01 

7.9236645E-04 

5 

5 

1 .9999999E-01 

1. 1200621E-03 

5 

6 

2.5000000E-01 

1.4795458E-03 

11 

7 

2.9999995E-01 

I.8708138E-03 

5 

8 

3.4999996E-01 

2.2938703E-03 

5 

9 

3.9999998E-01 

3. 7267059E-03 

5 

10 

4.4999999E-01 

3.8873379E-03 

5 

11 

5 .OOOOOOOE-01 

3. 7537606E-03 

11 

12 

5.4999995E-01 

5.907221  9E-05 

17 

13 

5.9999996E-01 

7. 759B752E-04 

11 

14 

6.4999998E-01 

1.0814792E-03 

12 

15 

6.9999999E-01 

2. 7733436E-04 

17 

16 

7.5000000E-01 

1.9064C1.  3E-03 

15 

17 

7.9999995E-01 

3.3322402E-05 

15 

18 

8.4999996E-01 

6.8813126E-05 

21 

19 

3.9999998E-01 

1.5895028E-05 

32 

20 

9.4999999E-01 

2.01  78050E-01 

5 

21 

l.OOOOOOOE  00 

6.0732605E-01 

11 

22 

1.04999926  00 

1.2206039E  00 

75 

23 

1.0999994E  00 

1.227S671E  00 

40 

24 

1.1499996E  00 

1.2347345E  00 

5 

25 

1.1999998E  00 

1.2412090E  00 

5 

26 

1.2500000E  00 

1.2472906E  00 

19 

27 

1.2999992E  00 

1.2535877E  00 

5 

28 

1.3499994E  00 

1.2597990E  00 

5 

29 

1.3999996E  00 

1.2659225E  00 

5 

30 

1.4499998E  00 

1.27163226  00 

5 

31 

1,50000006  00 

1.2779064E  00 

18 

32 

1.5499992E  00 

1.23425316  00 

5 

33 

1.5999994E  00 

1.2906771E  00 

5 

34 

1.6499996E  00 

1.2981052E  00 

5 

35 

1.699999BE  00 

1.3043680E  00 

5 

36 

1.7500000E  00 

1.31039816  00 

22 

37 

1.7999992E  00 

1,31619076  00 

5 

38 

1.B499994E  .00 

1.3217497E  00 

5 

39 

1.8999996E  00 

1.32659056  00 

5 

43 

1.9499998E  00 

1.3318396E  00 

5 

41 

2.0000000E  00 

1.3370161E  00 

11 

42 

2.0499992E  00 

1.34211646  00 

5 

Figure  4.  Partial  output  of  example 

2. 
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Because  of  the  time  required  to  compute  the  needed  Bessel 
functions,  the  call  to  POINTS  will  specify  101  points,  10  of  which  are 
initially  found  by  integration  and  one  point,  F(0)  = 0,  is  supplied.  The 
main  program  is  otherwise  identical  to  example  2.  IPRINT  = 1 was  used 
to  print  the  points  as  they  were  computed. 

The  transform  routine  is  written  to  reflect  equations  (2), 
(32),  and  (33): 

SUBROUTINE  TRANS (S,Z) 

COMPLEX  S,Z,F,G,E,J0,J1 
DATA  S71/. 84261498/ 

F = (0.,1.)*CSQRT(8.*S) 

CALL  CBES01(F,J0,J1) 

E = l./(l.  - 2.*J1/(F*J0)) 

F = F*S71 

CALL  CBES01(F,J0,J1) 

22 


G = 1.  + .8*J1/(F*J0) 
G = S*CSQRT(E*G) 

E = CEXP(G) 

Z = 2./(S*(E  + l./E)) 

RETURN 

END 


) 


where  siibroutine  CBESOl (F, JO, Jl)  was  a local  program  to  compute  the 
Bessel  functions  Jq(F)  and  Ji (F)  by  a series  expansion. 

Figure  6 shows  a partial  output  of  the  computed  points,  taking 
just  24  s of  CPU  time  on  an  IBM  360/195.  All  the  points  not  shown  were 
computed  by  interpolation. 


T IKF  FUNCTION  tRRCR  CODE 


1 

0.0 

0 .0 

1 

11 

l.OOOOCOO? 

00 

-t . 1 16E224E-03 

1 1 

21 

2.3300000E 

33 

1.1911774E  00 

12 

31 

3.0000000E 

00 

1.3226652E  00 

16 

41 

4 .OOOOOOOE 

00 

1.0734C15E  00 

12 

51 

5 .OOOOCOOE 

00 

9.4529492C-01 

14 

61 

6,  OOOOOOOE 

00 

9.368028CE-01 

17 

71 

7. OOOOOOOE 

00 

9.B935735E-01 

16 

81 

8. OOOOOOOE 

00 

1 .0182199E  00 

18 

91 

9, OOOOOOOE 

00 

1.012173  7E  00 

21 

101 

1 .OOOOOOOE 

01 

1.00015640  00 

23 

6 

5.0000000E- 

•01 

5.9992250E-05 

15 

16 

1.5C00003E 

00 

1.319S050E  00 

11 

26 

2.50C00OCE 

00 

1.2702408E  00 

11 

36 

3.5303333E 

33 

1.2153699E  33 

14 

56 

5.5300030E 

00 

9.1 71983CE-01 

16 

66 

6.5000000E 

00 

9.655851 lE-Cl 

13 

IS 

1 .6999993E 

33 

1.1131294E  00 

12 

23 

2.1 999998F 

00 

1.2279797E  03 

11 

28 

2.5999998E 

00 

1.2932730E  00 

14 

33 

3. 1999998E 

00 

1.322&388E  OC 

16 

33 

3 .6999993E 

0 3 

1.14847566  00 

15 

43 

4.1999998E 

03 

1.0367432E  03 

12 

53 

5. 199999SE 

00 

9.2949528E-01 

1 5 

9 

7.99°9995E- 

•01 

3.9430073E-05 

11 

12 

1.3999994E 

33 

4,21337016-01 

32 

14 

1.2999992E 

00 

8.55967586-01 

14 

17 

1 .59999946 

00 

1.0704C69E  00 

11 

1 3 

9.9999998E- 

31 

3. 5557676E-35 

22 

2 

9.9999964E- 

32 

2.  14246S3E-05 

5 

3 

1.9999999E- 

■01 

3.81 36262E-05 

5 

4 

2.9999995E- 

•31 

5.3134695E-35 

5 

5 

3.9999998E- 

•01 

5. 7419995E-05 

5 

7 

5 .9999996E- 

01 

5. 164735  2E-05 

5 

8 

6.9999999E- 

01 

4.4793269E-05 

5 

13 

l . 1 999998 E 

00 

7.0375639E-01 

19 

15 

1.3999996E 

00 

9.5463C85E-01 

5 

Figure  6.  Partial  output  of  example  3 (cont'd). 
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19 

1 .7999992E 

00 

1 . 1420450E 

CO 

5 

20 

1 .3999996,6 

00 

1.1685181E 

00 

5 

22 

2.0999994E 

00 

1.21086886 

00 

5 

2<i 

2.2999992E 

00 

1.2434921E 

03 

5 

25 

2.399Q996E 

00 

1.25747376 

00 

5 

27 

2. 59Q0994E 

03 

1.23227346 

33 

5 

29 

2.7999992F 

03 

1.3337567F 

33 

5 

30 

2 .8999996E 

00 

1.31455616 

03 

5 

32 

3.0999994F 

00 

1,32661826 

00 

5 

34 

3.2999992E 

00 

1 . 30176936 

00 

5 

35 

3. 3999996E 

33 

1.25374236 

33 

5 

37 

3 .5999994b 

03 

1 . 1 Bl 77996 

00 

5 

39 

3 .7999992F 

00 

1. 12031S4E 

CO 

5 

40 

3.3999996E 

33 

1.39543346 

00 

5 

Figure  6.  Partial  output  of  exeunple  3. 


Figure  7 shows  a graph  of  the  results.  The  smoothness  of  the 
computed  points  suggests  that  good  graphical  accuracy  has  been  obtained. 
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1.3  - 
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0 


L 
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TIME  (S) 


Figure  7.  Graph  of  example  3. 
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6.  C0NCLUSI.*5NS 


On  the  basis  of  the  three  examples  given  here  and  others  tried, 
subroutine  POINTS,  in  conjunction  with  others  listed  in  the  appendix,  is 
am  effective  method  for  inverting  a Laplace  transform.  The  adaptive 
interpolative  procedure  helps  cut  computer  time  by  not  requiring  all 
points  to  be  computed  by  directly  computing  the  Bromwich  integral. 
Thus,  costs  are  nominal  even  for  relatively  complicated  tramsforms. 
Since  the  use  of  the  program  is  relatively  easy,  it  is  deemed  a useful 
asset  to  any  software  library. 
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APPENDIX  A.— PROGRAM  LISTING 


This  appendix  contains  the  listing  of  the  subroutine  POINTS  and  all 
subroutines  it  calls  with  the  exception  of  TRANS.  The  user  must  write  a 
main  program  that  calls  POINTS  (or  TPOINT  if  only  one  point  is  to  be 
calculated)  and  the  subroutine  TRANS. 
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SUORCDf INe  P0l^rs( TRANS, P,  F.M.Wl,  N I , A? , H , V 1 , t P , I PR  I M I 

rms  KduriNt  coKPurrs  rut  iNvtRsc  laplace  transecrm  fcr  a sei  cf 
pcnns  specified  by  ehe  user,  the  ARCUPENr  list  is  ocfinco  as 

F-CLLOWS  . 

TRANS  - THE  POP.'-IY  NAME  OF  THE  SUBROUTINE  THAT  COMPUTES  THE 

TRANSFER  FUNCTION  VALUES,  AND  ViHICH  MUSI  RE  WRIIIEN  BY 
THE  USER.  ITS  EXACT  FORM  IS 
SUBROUTINE  TRANSIS.FI 
CCMPLEX  S,F 

' FOR  ANY  COMPLEX  INPUT  S,  F = FISI  MUST  BE  CALCULATED 

WITHOUT  DESfROTING  S.  THE  ACTUAL  NAME  CF  THIS  SUBROUTINE 
MUST  CE  DECLARED  TO  BE  EXTERNAL  BY  THE  CALLING  PROGRAM. 

P - THE  MAXIMUM  CF  THE  REAL  PARTS  OF  ALL  SINGULARITIES  OF 

TRANS.  THESE  INCLUDE  POLES,  BRANCH  POINTS,  AND  ESSENTIAL 
SINGULARIT  lES. 

E - THE  ABSOLUTE  ERROR  DESIRED  FOR  THE  POINTS  CALCULATED. 

THIS  IS  MORE  OF  A GUIDE,  SINCE  THE  ACTUAL  ERROR  MAY  BE 
SI  I GHTLY  GREATER . 

M - THE  MAXIMUM  NUMBER  OF  2.^^P1  INTERVALS  CCNSICERED.  IF  FIT! 
IS  CONTINUOUS  WITH  CONTINUOUS  DERIVATIVES,  133  IS  USUALLY 
SUFL  IC  lENT  . 

Wl  - A WORKSPACE  CIM.ENSIONED  M IN  THE  CALLING  PROGRAM. 

Nl  - THE  DIMENSION  OF  Tl,  Yl,  AND  IE  BELCH.  IMS  IS  THE  NUMBER 
OF  POINTS  TO  BE  COMPUTED.  SHOULD  BE  .GE.2  AND  .LE.501 
N2  - THE  NUMBER  OF  POINTS  INITIALLY  COMPUIEC  BY  INTEGRATION 
OF  THE  BROMWICH  INTEGRAL  IN2.GE.2J.  IF  THE  FUNCIICN  IS 
SMOOTH,  N2  = 5 OR  I 0 IS  USUALLY  SUFF IC IE  NT . THE  PCTNTS 
COMPUTED  ARE  APPROXIMATELY  ECU! SPACED  OVER  Nl.’ 

Tl  - THE  ARRAY  OF  TIMES  FOR  WHICH  THE  INVERSE  TRANSFORM  IS 

DESIRED.  THE  VALUES  OF  Tl  MUST  BE  STORED  PRIOR  TO  CALLING 
POINTS  AND  ORDERED  TO  INCREASE  WITH  I. 

Yl  - THE  ARRAY  CF  ANSWERS.  Yllll  CORRESPONDS  TO  Till). 

IE  - INPUT/OUTPUT  ARRAY.  AS  INPUT, 

lEIIl  = C MEANS  THAT  Yl(I)  MUST  BE  COMPUTED  BY  THE 
PROGRAM  AS  IT  SEES  FIT. 

= I MEANS  THAT  Yllll  IS  SUPPLIED.  ClF  Till) 

• LE.O,  Yim  MUST  BE  SUPPLIED. I 
= 3 IS  A REQUEST  TO  COMPUTE  Yllll  VIA  THE 
BROMWICH  INTEGRAL. 

AS  OUTPUT, 

lEII)  = L I.GE.Il)  MEANS  VKII  WAS  SUCCESSFULLY 

COMPUTED  VIA  THE  BRUMWICH  INTEGRAL  WITH  L 
INTERVALS  UF  LENGTH  2.«PI. 

= 5 MEANS  Ylin  WAS  COMPUTED  BY  INTERPOLATION. 

= 2 MEANS  Ylin  HAS  NOT  FOUND  TO  THE  ACCURACY 
DESIRED  IN  M INTERVALS  OF  LENGTH  2.»P1.  THE 
BEST  POSSIBLE  ANSWER  IS  RETURNED. 

= I MEANS  THAT  Yllll  HAS  SUPPLIED. 
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= -L  Ml ANS  niAI  THE  COMPUIAflfN  OF  POINT  1 WAS 
ACORfEO  HECAOSE  THE  t TM  IMERVA.L  CCOIO  NOT 
UE  INrEGRAlEO 

IPRINI  - THE  PRINT  CONTROL. 

IF  IPRINT  = 0,  MO  OUTPUT  IS  PRINTED. 

= I,  PRINTED  AS  POINTS  ARE  CCPPUTOO  INOT 
ORDERED! . 

= 2,  PRINTING  AFTER  ALL  PCINIS  ARE  COMPUTED 
lOROEREOI . 

ANV  POINT  COMPUTED  BY  INTEGRATING  THE  BROMWICH  INTEGRAL  IS  DONE  BY 
CALLING  TPOINT. 

INTERPOLATION  IS  DONE  BY  SUtROUTIME  INTERP,  WHICH  IS  A RCDIFIEO 
ROUTINE  FROM  THE  IBM  360  SSP. 


DIMENSION  Willi,  Til  II 

* VALI501) 

EXTERNALTRANS 
PRINT  HEADER  IF  REOUIPEC. 

lFIIPRlNT.GE.il  WR1TET6,22> 


ni  II,  Yll  11  , lEI  1 1 ,T(  SOI  I , YISOl  I ,ARGI  SOI  I , 


SET  UP  THE  Nl  POINTS  TO  BE  CALCULATED  BY  INTEGRATION.  J IS  THE  INDEX 
OF  THE  POINT,  FIRST  COMPUTED  IN  FLOATING  POINT  BY  FJ,  AND  INCREMENTED 
GY  XN. 

F J = 1 . 

J=l 

<N=FLOAT  (Nl-n  / FLOAT  (N2-1  I 
3 DO  S 1=1, N2 

JFIT.fC.N2l  J=N1 
1F(  IElJl.EO.il  GOTO  7 
IFdETJl.NE.O)  GOTO  6 
IF  IT  IIJ  I .GT .0. I GOTO  ^ 
lEI  Jl=-l 
GOTO  7 

A CALL  TPOINT (TRANS, P,T1IJI . Y 1 ( J I , E . M, I E ( J I , k 1 1 

IFMEIJI.EO.OI  IE(JI  = 2 

7 IFIlPRINT.EO.il  WR1TEI6.2A)  J , T 1 1 J I , Y 1 (J I , I E I J I 
6 FJ=FJ»XN 

J=F J».S 
S CONTINUE 

COMPUTE  SPECIAL  POINTS  FEOUESIED  IlETJI  = 21. 

00  a J=1,N1 

IF  I lEl  J I .NE .31  GOTO  8 

IFITK  JI.LE.O.  I GOTO  8 

CALL  TPC1NT(TRANS,P,T1UI,  Y I (J  I , E . M , I E IJ  I , Ml  I 
IFI  lEI  J I.EO.OI  IEUi  = 2 

IF (1PR1NT.EQ.1I  WR11Elb,2AI  J , T 1 1 J | , Y 1 1 J I , I E (J I 

8 CONI INUE 

FIND  INDICES  II  AND  12  SURROUNDING  UNCOMPUTED  POINTS. 
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9 12=1 
NN  = l 
13  J=3 

DO  20  1 = 12,  Nl 
IFU.NE.OI  GOTO  15 

IF(  lEI  n.NE  .O.ANO.(E(  .EQ.Ot  GOTO  18 
GO  TO  20 

15  IF(  (E(t  I.NE.O.  ANO.  lEI  I-n  .EQ.Ot  GOTO  17 
GOTO  23 

17  12=1 
GOTO  25 

18  J=l 
11  = 1 

20  CONTINUE 

if:i2.ne.i)  goto  9 

GOTO  70 

C CHOOSE  K ABOUT  MIDWAY  BETWEEN  II  ANO  12. 

25  K=ll  IM  21/2 

C 'lOVe  ALL  GOOD  POINTS  TO  T ANO  7 ANO  TRY  TO  INTERPOLATE  POINT  K. 
J = 0 

00  33  I=1,N1 

IF! 1 lElI  l.NE.li  .ANO.  (1E11I.LT.  91)  GOTO  30 
J=J*1 

T1JI=T1III 
Y1 J)=YIII) 

30  CONTINUE 

C J IS  THE  NUMBER  OF  POINTS.  IF  LESS  THAN  5,  CALL  TPOINT. 
IF1J.LT.5)  GOTO  53 

CALL  INTERP  IJ, T, Y, ri 1K»,YIIK>,E,IE IK) ,NN,ABG,VAL1 
IFI  lEIK  I.NE .01  GOTO  50 
1E1K1=5 
GOTO  10 

C POINT  COULO  NOT  BE  INTERPOLATED.  CALL  TPOINT. 

53  CALL  TPCINT  ITRANSfP.TlIKl,  Y II K ) ,E . M , I E I Kl ,K1 1 

IFI  lEIKl.EO.OI  IE(KI  = 2 

GO  IF  I IPRINT.EO. 1)  WR1TEI6,2A)  K , T 1 1 K ) , Y 1 ( Kt  , I E I K ) 

GOTO  10 

C LAST  PASS.  MOVE  ALL  GGCO  POINTS  TO  T ANO  t. 

70  J=0 

00  83  1=1. Nl 

IF  I 1 E 1 1 I .NE  .1 . ANO.  lEl I I .L T.6)  GOTO  80 
J=  J*l 

T( JI=T11I) 

YIJ)=Y1  1 1) 

80  CONTINUE 
NN=I 

00  90  K=l,Nl 

IF  I lEIK I .EQ. l.CR. lEIK) .GT .5)  GOTO  90 
IFIJ.LT.5)  GOTO  95 

CALL  INTERPIJ.T, Y.TllKl, VTEST.E.IERR.NN, ARG.VALI 
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1F(  IE<)R  .NE.OI  CIjTO  <)S 

Y1  IKMYIEST 

IEIK>s5 

GOTO  98  ...... 

95  1FIIE(KI.E0.2.CR.IE(K).LT.0I  GOTO  90 

CALL  TPOINTITRANS.P.TI  IK)  > Y 1 1 K ) , E , H , I E I K ) , Vi  U 

IFIIE(K).EQ.O)  IEIK)=2 

98  IFIIPRINT.EQ.II  WR11E(6,2<i)  K , T 1 1 K ) i Y I ( K I , 1 E ( K I 
90  CONTINUE 

IFI IPRINT.LE.l)  RETURN 

C PRINT  OUTPUT  IF  IPR INT  IS  GREATER  THAN  1. 

22  FORMAT  lull, 9X,AHTIMEi8Xi8HHUNCTIONi2X,10HERRCR  CCOE  /) 
00  23  1=1, Nl 

2 3 NRITEI6,2AI  I , T 1 1 I I , Y 1 1 I ) , IE  I I ) 

24  FORMAT  114, IP2E15.T, 15) 

RETURN 

END 


SUBROUTINE  TPC 1 N T I TR AN S ,P , T , Y ,EPS ,M , IE ,Wl ) 

TPOINT  COMPUTES  CNE  POINT  YIT)  OF  THE  LAPLACE  TRANSFORM  INVERSE  BY 
INTEGRATING  THE  BROMWICH  INTEGRAL.  THE  ARGUMENT  LIST  IS  AS  FOLLOWS. 

TRANS  - SEE  POINTS  FOR  DEFINITION 
P - SEE  POINTS  FOR  DEFINITION 

T - THE  TIME  FOR  WHICH  THE  INVERSE  LAPLACE  TRANSFORM  IS 
DESIRED 

Y - THE  ANSWER 

CPS  - THE  ABSOLUTE  ERROR  DESIRED.  THIS  IS  MORE  OF  A GUIDE, 

SINCE  THE  ACTUAL  ERROR  KAY  BE  SLIGHTL.  GREATER. 

M - SEE  POINTS  FCR  DEFINITION 
IE  - OUTPUT  ERROR  CODE. 

IE  = 0 MEANS  THAT  THE  ACCURACY  COULD  NCT  BE  OBTAINED 
IN  M INTERVALS  OF  WIDTH  2.*PI.  Y IS  THE  BEST 
ANSWER  OBTAINED. 

= L I.GT.  3)  MEANS  THAT  THE  ANSWER  WAS  OBTAINED 
IN  JUST  L 2.»PI  INTERVALS. 

= -L  l.LT.  3)  MEANS  THAT  THE  RUN  WAS  ABORTED 

BECAUSE  OF  DIFFICULTY  INTEGRATING  INTERVAL  L. 
W1  - WORKSPACE  OF  DIMENSION  M. 

INTEGRATION  OF  THE  BROMWICH  INTEGRAL  IS  DONE  BY  SUBROUTINE  SUM.  IT 
SUMS  THE  RESULTS  CF  INTEGRALS  FOR  INTERVALS  OF  2.*PI  AND  ATTEMPTS  TO 
ACCELERATE  THE  SUM  BY  CERTAIN  TRANSFORMATIONS.  THE  BROMWICH  INTEGRAL 
IS  COMPUTED  ALONG  A PATH  S = GC /T  * P WIHERE  GC  IS  INITIALLY  .31. 

THIS  IS  CLOSE  TO  A SINGULARITY  AND  CONVERGES  IN  SUM  RAPIDLY.  IF 
DIFFICULTY  IS  ENCOUNTERED  IN  ANY  INTERVAL  BECAUSE  GC  IS  TOO  SMALL, 

GC  IS  DOUBLED,  FOLLOWED  BY  ANOTHER  ATTEMPT  TO  COMPUTE  THE  INTEGRAL. 
EXTERNAL  TRANS 
COMMQN/TPT/G,Nl,E,T 1 
COMMQN/O/PI ,PI 2, XX, YY 
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OIHENSICN  Mini 
IHP.LT.II  KfcTUKN 
Pl=?.m593 
Pl2=6.a83l8S 
1E  = -1 

n = r 

GC  = . 0 1 

C Nl  CONIKOLS  THE  GAUSSUN  FOPKULA  USED  FOR  FHE  FIRST  INTEKVAl.  SET  TO 
C 5 TO  USE  A 32-POINT  FORMULA. 

Nl=5 

2 G=GC/I*P 

C ABORT  IF  G*T  IS  TOO  LARGE. 

IF  (G«T.GT.Il.)  return 

C EGT  IS  THE  COMMON  FACTOR  OF  ALL  TERNS  SUMMED. 

EGT=2.#EXPIG*T I/T 

C ADJUST  THE  ERROR  REQUIREMENT  FOR  EGT  AND  PASS  E IN  COMMON. 

E=£PS/EGT 

CALL  SUM!  V,M,E,  IE, TRANS, K.WII 

IFIIE.LT.O)  GOTO  A 
IFIIE.GT.OI  GOTO  3 

C ANSWER  IS  CK,  SET  IE  TO  K,  THE  NUMBER  OF  INTERVALS  NEEDED. 

IE=K 

I Y=Y*EGT 
RETURN 

C ANSWER  NOT  FOUND  WITHIN  M INTERVALS.  FLAG  WITH  IE=0. 

3 IE=0 
GOTO  1 

C SUM  COULD  NOT  INTEGRATE  ONE  OF  THE  INTERVALS.  DOUBLE  CG  AND  TRY  AGAIN 

4 GC=GC»GC 
GOTO  2 
END 


SUBROUTINE  SUM  I Y ,M , E , 1 E , TR ANS , K , X I 

SUM  COMPUTES  THE  BROMWICH  INTEGRAL  BY  OBTAINING  A SEQUENCE  OF 
PARTIAL  SUMS  AND  APPLYING  THE  DELIA  SQUARE  TR ANSF CRPA T I CN 

1.  TO  THE  SEQUENCE,  IF  IT  MONOICNIC  WHERE  THE  MAGNITUDE  OF  THE 
DIFFERENCE  IS  DECREASING  10  CONSECUTIVE  TIMES 

2.  TO  A SUBSEQUENCE  OF  PEAKS  WHERE  THE  HAXIMUMS  ARE  DECREASING 
AND  MIMMUMS  ARE  INCREASING 

3.  TO  A SUBSEQUENCE  OF  AN  ENVELOPE  OF  THE  PEAKS  WHERE  THE 
HAXIMUMS  ARE  INCREASING  OR  THE  MIMMUMS  ARE  DECREASING 

IF  3 CONSECUTIVE  PROJECTIONS  ARE  WITHIN  E OF  EACH  OTHER,  THE  SUM  IS 
CONSIDERED  FOUND.  IF  THE  PROJECTIONS  OSCILLATE,  A DELTA  SQUARE 
TRANSFORMATION  IS  APPLIED  TO  THEM  TO  OBTAIN  A FINAL  ANSWER.  OTHER- 
WISE THE  LAST  PROJECTION  IS  ACCEPTED  AS  THE  FINAL  ANSWER, 

THE  ARGUMENT  LIST  IS  AS  FOLLOWS 

Y - THE  BROMWICH  INTEGRAL  ANSWER  (WITHOUT  THE  FACTOR  EGT  IN 
TPOINTl  . 

M - SEE  POINTS  FOR  DEFINITION 
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f.  - SCE  POINTS  TOR  OtFINITION 
IE  - OUTPUT  ERROR  COOH 

le  = 0 lEANS  THAI  y WAS  FCUNO 

= M MEANS  THAT  TOE  ACCURACT  COULD  NOT  BE  ACHIEVED.  THE 
LAST  VALUE  OF  THE  PARTIAL  SUM,  XIH),  IS  RETURNED 
AS  THE  final  answer. 

= -L  l.LT.  01  MEANS  THAT  THE  RUN  WAS  ABCRTEO  BECAUSE 
OF  DIFFICULTIES  IN  INTEGRATING  INTERVAL  L. 

TRANS  - SEE  POINTS  FCR  DEFINITION 

K - THE  NUMBER  OF  INTERVALS  USED  IF  IE  = 0. 

X - THE  WORKSPACE  OF  DIMENSION  M,  USED  TO  STORE  THE  PARTIAL 
SUMS  OF  THE  BROMWICH  INTEGRAL. 

DIMENSION  XI  l),DEI3»,Plll»,Rim 
external  TRANS 
C < IS  THE  ARRAY  OF  PARTIAL  SUMS 

C NR  IS  THE  NUMBER  OF  ELEMENTS  IN  THE  R ARRAY  THAT  STORES  ENVELOPE 
C INFORMATION  OF  THE  PEAKS  CF  X. 

C NP  IS  THE  NUMBER  CF  ELEMENTS  IN  THE  P ARRAY  THAT  STORES  THE  PEAKS 
C OF  X. 

C NO  IS  THE  NUMBER  CF  CONSECUIIVE  VALUES  OF  X THAT  ARE  PCNOTOMC  WITH 
C DECREASING  DIFFERENCES 

NR=0 
NP=3 
N0=0 

C COMPUTE  THE  INTEGRALS  OF  FACH  INTERVAL  IN  A DO  LOOP. 

DO  105  l=l,H 
K=  I 

C 0 IS  THE  VALUE  OF  THE  INTEGRAL  OF  THE  CURRENT  INTERVAL. 

C 01  IS  THE  VALUE  OF  THE  PREVIOUS  INTERVAL. 

C DA  ANC  Dl  ARE  THE  ABSOLUTE  VALUES  OF  D AND  01. 

0=FCTI I, TRANS, IE) 

D2  = D 

C ABORT  IF  IHERE  IS  TROUBLE  IN  COMPUTING  0. 

IFIIE.lt. 0)  RETURN 
DA=ABSIDI 
lFtI.NE.il  GOTO  5 
XI 11=0 
GOTO  100 

C SET  UP  X VALUES  AND  COMPARE  DIFFERENCES. 

5 xin=xii-n  ♦ D 

IFIOA.GT.OLI  GOTO  10 

C THE  DIFFERENCE  IS  LESS  THAN  THE  PREVIOUS  DIFFERENCE, 

ND=NDM 

IFII.LT.llI  GOTO  too 
C CHECK  FOR  OSCILLATING  SEQUENCE. 

■IFt  0*01. LT.O.)  GOTO  10 
1F1ND.lt. 10)  GOTO  20 

C THERE  ARE  13  CONSECUTIVE  HONOTONIC  VALUES  OF  XIII.  EXTRAPOLATE  THE 
C LAST  5 VALUES  WITH  THE  DELTA  SQUARE  TP ANSFORMA T ICN. 
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CML  otLrAZ(5,  ,oe) 

C RtJCCr  IF  ALL  3 LXrRAPnLMIlO  VALUES  ARE  NOT  WITHIN  E CF  EACH  OTHER. 
IF  I ABS  lOEl  I l-OEUn  .GT  .E  I GOTO  133 
IF|ABS(CE(3I-0EI2I) .GT  .Ei  GOTO  100 
IFIA8S(0E(3I-0E(IM.GT.E)  GOTO  lOO 
11  1F(I0E13)-DEI2)»*IDE121-CE1U1  .gT.O.)  GOTO  8 
C EXTRAPOLATE  THE  EXTRAPOLATED  VALUES  IF  THEY  OSCILLATE. 

CALL  CELT A2 13. CEi CEI 

Y=DElll 

GO  TO  9 

C ANSWER  ACCEPTED. 

8 Y=0EI3I 

9 IE=0 
RETURN 

C THE  MAGNITUDE  OF  THE  DIFFERENCES  OF  THE  XII)  ARE  INCREASING.  CHECK 
C FCR  PEAKS 
10  N0=0 

IFII. LT. 11)  GOTO  100 
20  11=2 

IFINP.NE.O)  11=1PVI 
IMl=l-l 

00  15  J=1I, IHl 

C 0 AND  01  CHANGE  MEANINGS  ANO  ARE  NOW  TEMPORARY  VARIABLES. 

0=XI Jl-Xl J-l) 

01=XIJ*1)-<1J) 

IFIO.GE.O.)  GOTO  17 
1FI01.lt. 3.)  GOTO  15 

C A MINIMUM  IS  DETECTED.  REFINE  WITH  SUBROUTINE  PEAK. 

NP=NPH 
IP  = J 

CALL  PEAKIPINP),  <U)tO,Dl) 
lA  IFINP.LE.2I  GOTO  IS 

C REJECT  IF  MINIMUM  IS  LESS  THAN  THE  PREVIOUS  MINIMUM. 

IFIPINPI.LE.PlNP-2))  GOTO  16 
C REJECT  IF  XIJ*1)  IS  .GE.  THE  PREVIOUS  MINIMUM. 

IFIXIJ*1I.LT.PINP-1I)  GO  TO  18 
GOTO  23 

C THE  PEAK  IS  REJECTED.  STORE  IN  THE  R ARRAY. 

16  NP=NP-2 
NK=NR*1 

RINR)=|P|NP)*PINPUI  )*.5 

1 J=1 

PINP)=P(NP*2) 

IFINR.EC.II)  GOTO  25 

C CONTINUE  CHECKING  THE  NEW  PEAK  AGAINST  THE  PREVIOUS  ONES. 

GOTO  14 

18  IFINP.lt. 5)  GOTO  15 

IFIJ. NE.IMl)  GOTO  23 

C THERE  ARE  AT  LEAST  5 PEAKS.  COMPUTE  ANO  CHECK  THE  DELTA  SQUARE 
C TRANSFORMATION. 
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CALI  CrLTA2l5,PlNP-<i),OE) 

C CHECK  HIE  DlfFEKEKCE  OF  f RAN SFURMEO  VALUES. 

iriABSICEI  n-OE(2l  ».LE.E»  GOTO  21 
C ERRORS  TOO  GREAT.  .'OVE  PEAKS  IN  P ARRAT  I F NP  IS  .CE  . ll. 

23  IFINP.LE.IOI  GCrO  15 
DO  22  L=l. 10 
22  PILJ=PIL*U 
NP  = 10 
GOTO  15 

C CONTINUE  CHECKING  THE  DIFFERENCES  OF  THE  TRANSFORMED  VALUES. 

21  IF(AdS(CEI2i-DE(3l).GT.E)  GOTO  23 
IFIAHSIDEm-DE13»).GT.F»  GOTO  23 
C THE  ERROR  CRITEREA  ARE  SATISFIED  SO  THE  ANSWER  IS  FCtNO. 

GOTO  ll 

17  IFIDI.GT.O.)  GOTO  15 

C A MAXIMUM  IS  DETECTED.  REFINE  WITH  SUBROUTINE  PEAK. 

NP  = NP»-l 
IP  = J 

CALL  PEAKIPINP),XIJ),0,D11 
19  IFINP.LE.2)  GOTO  15 

C REJECT  IF  THE  MAXIMUM  IS  GREATER  THAN  THE  PREVIOUS  ONE. 

IFIPINP).GE.P(NP-2I»  GOTO  2A 
C REJECT  IF  XlJfll  IS  .LE.  THE  PREVIOUS  MINIMUM. 
IFIXUfll.GT.PiNP-in  GOTO  18 
GOTO  23 

C THE  PEAK  IS  REJECTED.  STORE  IN  THE  R ARRAY. 

2A  NP=NP-2 
NR=NR*l 

RINRI=IPINPI*PINP+1I )*.5 
IJ=2 

PINPI=PINP*2) 
lFINR.EQ.il  I GOTO  25 

C CONTINUE  CHECKING  THE  NEW  MAXIMUM  AGAINST  THE  PREVIOUS  ONES. 
GOIC  19 

C CHECK  R ARRAY.  IF  THE  MAGNITUDE  OF  THE  DIFFERENCES  IS  MCNOTONIC 
C DECREASING,  THEN  EXTRAPOLATE. 

25  DO  26  L=l,9 

IFIABSIR|Lf2)-R|L*^in.GT.ABSIRlL*l)-RILm  GOTO  27 

26  CONTINUE 

C THE  R ARRAY  VALUES  SEEM  TO  BE  APPROACHING  A LIMIT. 

CALL  0ELTA2I5,RI7I,0E) 

C CHECK  THE  DIFFERENCES  OF  THE  EXTRAPOLATED  VALUES  FOR  ERROR. 
IFiABSI0£m-CEI2i>.GE.EI  GOTO  27 
IFIABSIDEI3I-DEI 21  I.GE.E)  GOTO  27 
IFIABSI0EI3l-DEiin.GE.EI  GOTO  27 
C THE  ERROR  CRITEREA  ARE  SATISFIED  SO  THE  ANSWER  IS  ACCEPIEO. 

GOTO  II 

C R ARRAY  EXTRAPOLATION  YIELDS  TOO  MUCH  ERROR,  MAKE  ROOM  FOR  NEXT 
C VALUE. 

27  DO  28  L=l, 10 
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28  R(U 
NR  = 10 

CCNFINUE  SEARCHING  PCR  PEAKS  IN  LOOP. 

15  CONTINUE 

UPOME  AND  CONTINUE  MAJOR  LOUP. 

100  OL=OA 
10*5  01=02 

THE  LIMIT  Of  THE  NUN3ER  OF  INTERVALS  IS  REACHED.  RETURN  WITH  MIH 
PARTIAL  SUM. 

Y=X(M) 

IE=M 

RETURN 

END 

SUPROUTINE  INTERP 
PURPOSE 

TO  INTERPOLATE  FUNCTION  VALUE  V FOR  A GIVEN  ARGUMENT  VALUE 
X USING  A GIVEN  TABLE  IXARG.YVAL)  OF  ARGUMENT  AND  FUNCTION 
VALUES. 


USAGE 

CALL  INT  ERPlNOIN.XARGf  YVAL , X, Y.EPSi lER ,NN ,ARG .VAL» 


DESCRIPTION  OF  PARAMETERS 


AN  INPUT  value  which  SPECIFIES  THE  NUMBER  OF 
POINTS  IN  TABLE  IXARG.YVAL). 

THE  INPUT  VECTOR  lOIMENSION  NOIM.)  CF  ARGUMENT 
VALUES  OF  THE  TABLE  INOT  OESIRCYEOI. 

THE  INPUT  VECTOR  IOIMEN*^ION  NOIM)  OF  FUNCTION 
VALUES  OF  THE  TABLE  INOT  OESIROYED). 

THE  ARGUMENT  VALUE  SPECIFIED  BY  INPUT. 

THE  RESULTING  INTERPOLATED  FUNCTION  VALUE. 

AN  INPUT  CCNSTANT  WHICH  fS  USED  AS  UPPER  BOUND 
FOR  THE  ABSOLUTE  ERROR. 

A RESULTING  ERROR  PARAMETER. 

AN  INPUT  ESTIMATE  FOR  THE  SUBSCRIPT  TO  MAKE 
ABSIXARGINN)  - X)  A MINIMUM.  NN  RETURNS  SET  TO  THE 
CORRECT  VALUE. 

WORKSPACE  OF  DIMENSION  NDIH. 

WORKSPACE  OF  DIMENSION  NDIM. 


REMARKS  i 

ID  TABLE  IXARGfYVALi  SHOULD  REPRESENT  A SINGLE -V ALUEO 

FUNCTION  AND  SHOULD  BE  STORED  SO  THAT  XARGI I ) . GE . XARG IJ) 
WHET'WER  l.GT.J. 

12)  NO  ACTION  BESIDES  ERROR  MESSAGE  IN  CASE  NOIM  LESS 
THAN  1. 

131  INTERPOLATION  IS  TERMINATED  EITHER  IF  THE  DIFFERENCE 
BETWEEN  TWO  SUCCESSIVE  INTERPOLATED  VALUES  IS 
ABSOLUTELY  LESS  THAN  TOLERANCE  EPS.  OR  IF  THE  ABSOLUTE 
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VAIUE  OF  IMIS  r>l  rFCRFNCE  STOPS  0 1 1'.  IN  I SlU  NG  , OR  -VFICR  AL  I 

(NOIM-II  STEPS.  FURTHER  IT  IS  TERPIKMF.O  IF  THE  ALT 

PROCEDURE  DISCOVERS  rwO  ARGUMENT  VALUES  IN  VECICR  APG  AL I 
WHICH  ARE  IDENTICAL.  OEPENOENT  ON  THESE  FOUR  CASES.  AL I 

ERROR  parameter  IER  IS  CODED  IN  THE  FOLLCWING  FORM  ALT 

1ER=3  - IT  WAS  POSSIBLE  TO  BEACH  THE  BECUIREO  AL I 

ACCURACY  INO  ERROR!.  AL I 

I£R=l  - IT  WAS  I.1POSSIOLE  TO  REACH  THE  PECUIREO  AL  I 

ACCURACY  RECAUSE  OF  ROUNDING  ERRCRS.  INCREASE  AH 
EPS  ANO/OR  THE  ACCURACY  OF  THE  TABLE. 

I£R=2  - IT  WAS  IHPOSSIRLE  TO  CHECK  ACCURACY  BECAUSE  ALT 
NOIM  IS  LESS  THAN  3,  OR  THE  RECU I RED  ACCURACY  ALT 
COULD  NOT  BE  REACHED  BY  MEANS  OF  THE  GIVEN  ALT 

TABLE.  NDIH  SHOULD  BE  INCREASED. 

IER=3  - THE  PROCEDURE  DISCOVERED  TWO  ARGUMENT  VALUES  AL I 
IN  VECTOR  ARC  WHICH  <'RE  IDENTICAL.  ALT 

ALI 

SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED  ALI 

NONE  AL I 

ALI 

METHOD  ALT 

INTERPOLATION  IS  DONE  BY  MEANS  OF  AITKENS  SCHEME  OF  ALT 

LAGRANGE  INTERPOLATION,  ON  RETURN  Y CONTAINS  AN  I N I ERPOL AT EOAL I 
FUNCTION  VALUE  AT  POINT  X,  WHICH  IS  IN  THE  SENSE  CF  REMARK  ALI 
I3I  OPTIMAL  WITH  RESPECT  TO  GIVEN  TABLE.  FOR  REFERENCE,  SEE  ALI 
F.B.HILDEBRANO,  INTRODUCTION  TO  NUMERICAL  ANALYSIS.  ALI 

KCGRAW'HILL,  NEW  YORK/ TORONIO/LONDON . 1956,  PP.A9-50.  ALI 

ALI 

SU3ROUT  INE  INT ERP I NO IH , XARG, YVAL . X , Y . EPS . I ER ,NN , ARG , VAL I 


ALI 

DIMENSION  XARGl  1 1 . YVAL ( II , ARGI II , VALI 1) 

ALI 

IER»2  All 

IFINCIM-H9.37.31 

FIND  THE  CORRECT  VALUE  OF  NN 
I l-NN 

DELT2=X-XARGtlII 

START  SEARCH  LOOP 
00  25  K=1,NDIM 
JJ=II  ' 

IF(0eLr2I22,23,24 

IF  0ELT2.E0.0  THEN  NN  « II  AND  Y = YVALINN)  EXACTLY 
NN=  1 1 

Y=yVAlINNI 

1ER«0 

RETURN 
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IF  oti.r2.Gr.o  riiGN  rp.v  larger  suascRipr 

2«  I I = 1 I H. 

I:  ( I f .LE.NOIMI  GO  TO  26 
JJ  =N01M 
GO  TO  ZT 

IF  OCLTZ.LT.O  THEN  TRY  SRALLER  SOdSCRIPT 
ZZ  11=11-1 

IFdl.GE.n  GO  TO  Zb 
JJ  = l 

GO  TO  27 
Z6  OELTl=DELIZ 

OELTZ'X-XAPGII I ) 

COMPARE  DElIZ  KITH  DELTl.  IF  GREATERf  THEN  CLOSEST  POINT  IS  PAST 
IFIAESICELTZI.GE.ABSIOELTll)  GO  TO  27 
25  COMINUE 
27  NN=JJ 
H = JJ 

ARGI I l=XARG III) 

VAl I 1 )=yVAl I I I I 
OELTZ=0. 

J=1 

TRANSFER  TO  DETERMINE  THE  SECOND  CLOSEST  POINT. 

GO  TO  6 


START  OF  AITKEN-LOCP 

1 OELTl=DELT2 
I£NO=J-l 

00  Z 1 = 1.  TEND 
H=ARGI  I ) -APGU) 

IFIHIZ, 13,2 

2 VALIJ)  = IVALU1*(X-APG(JI  l-VAU  J I*  I X- ARGI  I 1)1 /H 
0ELT2  = ABSI,VAL1  J)  -VAL  I I END)  I 

IF  IJ-3)6,3,3 

3 IFI0ELT2-EPS)10,10,A 
A IF(J-5)6,5,5 

5 IFIDELT2-0ELT1)6, 11, 11 

6 J=J*1 

END  OF  AITKEN-LOOP  BUT  ME  MUST  FIND  THE  JTH  CLOSEST  POINT  BEFORE 
LOOPING  BACK  TO  STATEMENT  1. 

TFT J.GT. NOIN)  GO  TO  36 
IF(II.EO.l)  GO  TO  30 
IFIJJ.EC.NO)H)  GO  TO  29 

IFI  ABSIXARGI  Il-D-X)  .GT.  ABSI  XARCIJJ  *1)-X)  ) GO  TO  30 
29  11=11-1 

ARG(J)=XARGIII) 

VALI  J)  = yVALIll ) 


38 


oo  r>o  oo  on 


I'  W /III  • ■ 


CO  ro  i 

30  JJ=JJ»l 

ARGt J)=XARG(JJ I 
VALUJ=YVftL(  JJ) 

GO  ro  I 

DEFINE  VAL(U  IN  CASE  NOIH  = 1. 

37  VALt  1)  = YVALI U 
36  J=ND1M 

8 Y^VALIJI 

9 RETURN 

THERE  IS  SUFFICIENT  ACCURACY  WITHIN  NOIK-l  ITERATION  STEPS 

10  IER=0 
GOTO  8 

TEST  VALUE  0ELT2  STARTS  OSCILLATING 

11  IER=l 

12  J=IENO 
GOTO  8 

THERE  ARE  TWO  IDENTICAL  ARGUMENT  VALUES  IN  VECTOR  ARG 

13  1ER=3 
GOTO  12 
END 



FUNCTICN  Tlx, TRANS! 

C FUNCTION  T SETS  UP  THE  CORRECT  ARGUMENT  FOR  CALLING  TRANS  ANO 
C RETURNS  THE  PEAL  PART  OF  THE  COMPLEX  EVALUATION, 

COMPLEX  C 

COMMON/  TPT/G,Nl,E,n 

I CALL  TRANSICMPlXIG,X/ril,C) 

I=REAHCI 

RETURN 

END 



SUBROUTINE  PE AK I P, Y, D, D1 ) 

C PEAK  COMPUTES  THE  PEAK  OF  THE  3 POINTS,  Y-0,Y,  ANO  YtOl,  WHEN  IT  IS 
C KNOWN  THAT  0 AND  Dl  ARE  OF  OPPOSITE  SIGN.  A SIMPLE  2N0-CRDER 
C FORMULA  IS  USED. 

IF (Ol-D.NE.O. ) GOTO  5 
P=Y 

RETURN 

5 P=Y-I01*D1**2*. 125/(01-01 
RETURN 
END 



SUBROUTINE  OEL T A2 (N , X, Al 

C 0ELTA2  APPLIES  AITKEN'S  DELTA  SQUARED  TRANSFORMATION  TO  THE  N 
C ELEMENTS  OF  ARRAY  X.  THE  RESULTING  N-2  ELEMENTS  ARE  PLACED  IN  THE 
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AL  ( 

09  3 

ALI 

09A 

ALI 

095 

AL  I 

096 

ALI 

097 

ALI 

093 

ALI 

099 

ALI 

100 

ALI 

101 

ALI 

102 

ALI 

103 

ALI 

lOA 

ALI 

105 

ALI 

106 

ALI 

107 

ALI 

103 
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C 


c 


c 

c 

c 

c 

c 

c 

c 


c 


AP.i^AY 

OIKENSIDN  X(l)fA(l) 

NM2’N-2 

DO  5 1=1, NM2 

XI 1=X( I»ll 

XI2=XI  I*2>- 

o=xm-xii-xiuxi2 

IFI0.NE.3.)  GOIO  ^ 

A(  M = XI  I 
GOTO  5 

« Am  = XI2-(XI2-XI  n*»2/D 
5 CONTINUE 
RETURN 
END 


FUNCTION  FCTIN, TRANS, lERI 

FUNCTION  FCT  COMPUTES  TFE  BROMWICH  SUOINTERVAL  FRCM  2*PI*IN-U  TO 
2*Pl*N  BY  A HIGH-CRDER  GAUSSIAN  QUADRATURE.  THE  VALUE  IS  CHECKED  BY 
THE  NEXT  HIGHEST-CRJER  SUPPLIED.  PROVISION  IS  HADE  TO  USE  A LOWER 
PAIR  OF  QUADRATURE  FORMULAS  IF  THE  ERROR  IS  TOO  SMALL,  AND  A HIGHER 
PAIR  IF  THE  ERROR  IS  TOO  LARGE.  THE  ARGUMENT  LIST  IS  AS  FOLLOWS. 

N - THE  NUMBER  OF  THE  INTERVAL. 

TRANS  - POINTS  FOR  DEFINITION. 
lER  - ERROR  CODE 

lER  = 0 IF  FCT  IS  SATISFACTORY. 

= -N  IF  THE  ERROR  CRITERION  WAS  NOT  SATISFIED  WITH 
THE  30  AND  96  POINT  GAUSSIAN  QUADRATURE 
FORMULAS. 

EXTERNAL  TRANS 
COMMON/ TPT/G,N 1, EPS, Tl- 
COMMON/C/Pl ,PI 2, X,P 
DIMENSION  YI2I 

P IS  THE  MIDPOINT  OF  THE  INTERVAL. 

P=PI2*FL0AT INI -PI 

N2  CONTROLS  THE  ORDER  OF  THE  GAUSSIAN  FORMULA  SELECTED.  Nl  IS 
INITIALLY  SET  TO  5 IN  TPOINT 
N2=Nl 

K IS  THE  FLAG  TO  INDICATE  IF  THE  FIRST  OR  SECOND  FORMULA  IS  CHOSEN. 


K=l 

Nl»l=Il  MEANS  FAILURE 
AND  96-POlNT  GAUSSIAN 
WHERE  GC  as-  INCREASED 
13  NPl=Nl+l 

IFINP1.LT. Ill  GOTO 


TO  SATISFY  THE  ERROR  CRITERION  KITH  AN  80- 
QUADRATURE  FORMULA.  RETURN  IS  MADE  TO  TP''INT 
TO  GET  FURTHER  AWAY  FROM  SINGULARITIES. 

15 


IER=-N 


Nl=9 


RETURN 

SELECT  THE  APPROPRIATE  QUADRATURE  FORMULA. 
15  GOTOI 1 ,2,3.4,5,6, 7,8,9, 101 , N2 
1 CALL  Q6I  Y|KI,TRANSI 
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GOTO 

16 

2 

CALL 

0121 

Y IK  I , TRANS) 

GOTO 

16 

3 

CALL 

016  1 

Y(K),  TRANS) 

GOTO 

16 

A 

CALL 

02A  I 

YIK), TRANS) 

GOTO 

16 

5 

CALL 

0321 

YIK  I , TRANS) 

GOTO 

16 

6 

CALL 

QAOI 

YIK), TRANS) 

GOTO 

16 

7 

CALL 

0',8  1 

YIK), TRANS) 

GOTO 

16 

8 

CALL 

06  A I 

YIK),  TRANS) 

GOTO 

16 

9 

CALL 

0801 

Y IK  ),  TRANS) 

GOTO 

16 

10 

CALL 

096  ( 

YIK)  , TRANS) 

16 

IFIK 

.EO.  1) 

GOTO  20 

C CHECK  TO  SEE  IF  ER«OR  CRITERION  IS  SATISFIED. 

IFiaBSIvm-Yian.LE.EPS*.!  I GOTO  18 
C ERROR  CRITERION  IS  NOT  SATISFIED.  INCREASE  NI  AND  N2  AND  TRY  AGAIN. 
Y( l)=Y(2l 
Nl=Nl«-l 
N2  = NU1 
GOTO  13 

C ERRCR  CRITERION  IS  SATISFIED. 

18  FCT=YI2I 
IER  = 3 

C CHECK  TO  SEE  IF  ERROR  IS  TOO  SMALL. 

IF (ABSI VI  II-YI 21  ) .GT  .EPS*  I .E-31  RETURN  , 

C ERROR  IS  TOO  SMALL.  DECREASE  Nl  FOR  NEXT  INTERVAL. 

Nl=Nl-l 

IFIN1.lt. II  Nl  = l 
RETURN 
20  K=2 

N2=NPl 
GOTO  15 


SUBROUTINE  061  Y, TRANS) 

SUBROUTINE  06  COMPUTES  A-  6-POINT  GAUSSIAN  OUADRATURE  OVER  THE 
INTERVAL  AND  RETURNS  THE  ANSWER  Y.THE  ZEROS  OF  THE  LEGENDRE 
POLYNOMIAL  OF  DEGREE  6 FAVE  BEEN  MULTIPLIED  BY  PI  AND  STORED  IN  X, 
AND  THE  WEIGHTS  HAVE  BEEN  MULTIPLIED  BY  -COS  X AND  STORED  IN  U 
IBOTH  IN  DATA  STATEMENTS). 

COKMCN/C/Pl .PI  2. Z,P 
EXTERNAL  TRANS 
DIME-NSION  XI  3). Ml  3) 

DATA  X/2.929A39  ,2.077251 
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Y = 3. 

DO  5 I=li3 

2=xm 

5 y=Y*mPHtTRANS»  f HP-Z, raANS))»W(II 
RETURN 
ENO 



SUBROUTINE  0121  Y, TRANS) 

C SUBROUTINE  012  COMPUTES  A 12-POINT  GAUSSIAN  OUADRATURE  OVER  THE 
C INTERVAL  AND  RETURNS  THE  ANSWER  Y.  SEE  06  FOR  THE  MEANING  OF  X AND  H 
COMMCN/Q/PI ,PI2,Z,P 
EXTERNAL  TRANS 
DIMENSION  XI  6) ,WI  6) 

DATA  X/3.0B3664  ,2.343368  ,2.418721  ,1.843114  , 

♦ 1.155577  , .3934324  /,  W7 . 4 7C96 20  E- I , . 1 02 1 2 43  , 

♦ . 1200442  ,.5503602  E- 1 , - . 94 1 88 76  E - 1 , - . 2 301 1 1 9 / 

Y = 0. 

00  5 1=1,6 
Z=X( I ) 

5 Y=Yt(TlP*Z, TRANS)  ♦ T I P- Z , TR ANS ) ) I I ) 

RETURN 

ENO 



SUBROUTINE  0161  Y, TRANS) 

C SUBROUTINE  016  COMPUTES  A 16-POINT  GAUSSIAN  QUADRATURE  OVER  THE 
C INTERVAL  AND  RETURNS  THE  ANSWER  Y.  SEE  C6  FOR  THE  MEANING  OF  X AND  W 
EXTERNAL  TRANS 
C0HM0N/C/PI.PI2,Z,P 
DIMENSION  <1  8),W|  8) 

DATA  X/3. 108295  ,2.96747,2.719461  ,2.373173  , 

♦ 1.941115  ,1.433902  ,.8846836  ,.29849C6/,  W7. 2713741  E-l 

.6131218  E-l, .3683526  E - 1 , . 8963946E- 1 , . 54 1 437 3 E - I , -. 222461 3 E-l, 

♦ -.1 156855  ,-.1810734  / 

Y = 0. 

DO  5 1=1,8 

z=xm 

5 Y=YtlTIP«^Z, TRANS)  ♦ T ( P-Z , TR  ANS ))  *W  (I ) 

RETURN 

END 



SUBROUTINE  0241  Y.TPANS) 

C SUBROUTINE  024  COMPUTES  A 24-POINT  GAUSSIAN  QUADRATURE  OVER  THE 
C INTERVAL  AND  RETURNS  THE  ANSWER  V.  SEE  06  FOR  THE  MEANING  OF  X AND  H 
EXTERNAL  TRANS 
COMMON/0/PI,PI2,Z,P 
DIMENSION  XI12),HI12) 

DATA  X/3. 126473  ,3.062200  ,2.947676  ,2.784757  , 

*2.576112  ,2.325169  ,2.036046  ,1.713492  ,1.362802  ,.9897358  , 
*.6004176  ,.2012407  /,  W/1. 233982  E-2, 2. 844152  E-2, 
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£-^l  5.'i563l7i;-2,6. 192373  E- 2 , 5 . 902  5 7 3 E-2,^.379fc24  E-2i 
»1. 52  7986  E-2t -2.335162  E-2 , - 6 . 5 7062 3C - 2 1 - . I 038 265  ,-.  1253563  ! 

V = 0. 

00  5 1=1,12 
Z=X(I  ) 

5 Y=YMr(P»Z,rPANS»  ♦ TIP-Z,  TRANSn^wm 
RCTUFN 
END 

SUBRCUriNE  0321  Y,  TRAN'S! 

SUaflOUriNE  032  computes  a 32-point  GAUSSIAN  CUAOftATURE  OVER  THE 
INTERVAL  AND  RETURNS  THE  ANSWER  Y.  SEE  06  FOR  THE  MEANING  OF  X AND  W 
EXTERNAL  TRANS 
C07MGN/C/PI ,PI2,Z ,P 
DIMENSION  X I 16) ,WI 16) 

DATA  X73. 132997  .3.C96390  ,3,03089,2.937093  ,2.815876, 

♦2.668367  ,2.395934  ,2.300210,2.083015,1.836363  ,1,592373,1.323713 
*,1.032596, .7517333  ,.3536721  ,.1517630  /,  W/7. 018351  E- 

♦3,1.625777  E-2, 2.523663  E-2, 3. 355970  E-2,3.C58366  E-2, 3. 539352  E-2 

♦ ,3.587157  E-2,  3.385537  E-2, 3. 535757  E-2, 2. 127599  E-2, 1.805786  E-3, 
♦-2. 133759  E-2, -3.593979  E -2, -6 . 655 32 7 E-2 ,- 8 .5 555 69  E-2 , -9 . 53 3036 

♦ E-2/ 

Y = 0. 

DO  5 1=1,16 
Z = XI  I) 

5 Y=YMTIP*Z, TRANS)  ♦ T I P-Z  , TR  ANS ) ) *W  H ) 

RETURN 

END 

SUBROUTINE  C30I  Y, TRANS) 


SUBROUTINE  030 
INTERVAL  AND  R 
EXTERNAL  TR 
COKMCN/C/PI 
DIMENSION  X 
DATA 

♦ 2.530518  , 

♦ 2.111013  , 

♦ .3323239  , 
♦3,1.039383 
♦,3.539313  E 
♦2.122887  E- 
♦2,-3.973986 

Y = 0. 

DO  5 1=1,20 

z=xn  I 

5 Y=Y*mP*Z, 
RETURN 
END 


COMPUTES  A 30-POINT  GAUSSIAN  QUADRATURE  OVER  THE 
ETURNS  THE  ANSWER  Y.  SEE  06  FOR  THE  MEANING  OF  X AND  W 
ANS 

,PI2, Z,P 

(20),W( 20) 

X/3. 136056  ,3.112358  ,3.070153  ,3.009383  , 
2.833027, 2. 720392, 2. 59C596  ,2.335119  ,2.283938  , 
1.923395  ,1.726202  ,1.517627  ,1.299926  ,1.073306  , 
.6053773, .3636089  ,.1218071  /,  W/3. 521208  E- 

E-2, 1.637917  E-2, 2. 205171  E-2, 2. 731698  E-2, 3. 189002  E-2 
-2,3.737815  E- 2 , 3 . 7355 C IE- 2 , 3 . 386 25 8E-2 ,2 . 953 259  E-2, 
2,1.003032  F-2, -3.609110  E- 3 ,- 1 . 889 359  E-2 , -3 .371 256  E- 
E-2, -6. 258359  E- 2,-7 . 19 7 328  E -2 , - 7 . 693 168  E-2/ 


TRANS)  ♦ TI P-Z, TRANS) )*WII4 
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SUbPOUTINE  Y.rPANSJ 

C SUBRI1UIINE  0A8  CCMPUFCS  A Ae-POINf  GAUSSIAN  QUADRATURE  CVER  IME 
C INTERVAL  AND  RETURNS  THE  ANSWER  Y.  SEE  06  EOR  THE  R.CAMNG  CF  X AND  W 
EXTERNAL  TRANS 
CORMON/O/PI f PI 2i E,P 
DIMENSION  X12A|,WI2A) 

DATA  X/3.13T732  .3.121267  ,3.091710  .3.CA9203  , 

♦ 2.  993899  ,2.926338  ,2.3<t5903  ,2.  753832  , 2 .653211  ,2.535<V73  , 

♦ 2.AI0101  ,2.27',620  , 2.129599  ,l.9756A5  ,1.813905  ,1.693559  , 

♦ 1.956819  ,1.283926  ,1.055599  ,.9027759  ,.7061163  ,.5069950  , 

♦ .3097993  ,.1017253  /,  W/3. 153323  E-3, 7.326090  E-3, 

♦1.196296  E-2, 1.551237  E-2, 1.990260  E-2, 2. 302528  E-2 , 2. 62 3629E -2 , 

♦ 2.885332  E-2, 3.066299  E-2,  3. 192922  E-2, 3. 091698  E-2, 2. 391059  C-2, 

♦ 2.529998  E-2,  1.98359  1 E-2,  1.27361  1 E-2,9.C12932  E -3 , -5 .999829  E-3, 
«-l. 672663  E-2, -2. 7 7 7093  E- 2, -3 .892928  E-2 , -9  .302 28 1 E-2, 

♦ -5.589856  E-2, -6 . 1995 71  E-2, -6 .990303  E-2/ 

7=0. 

DO  5 1=1,29 
Z=XI  I 1 

5 V=Y»ITIP«-Z,TRANS1  ♦ T I P-7  , TR  ANSI  ) ♦»  I II 
RETURN 
ENO 

C 

SUBROUTINE  C69I  Y, TRANS! 

C SUBROUTINE  Q69  COMPUTES  A 69-POINT  GAUSSIAN  QUADRATURE  OVER  THE 
C INTERVAL  AND  RETURNS  THE  ANSWER  Y.  SEE  06  FOR  THE  MEANING  OF  X ANO  H 
EXTERNAL  TRANS 
CCMKON/Q/PI ,PI2,7,P 
DIMENSION  X(32),WI32! 

DATA  - - X/3. 139939  ,3.130095  ,3.113360  ,3.089292,- 

♦ 3.057796  , 3.019098  , 2 .97 3239, 2 .9 20 328  ,2.860990  ,2.793867  , 

♦ 2.723617  ,2.693915  ,2.559596  ,2.962922  ,2.365053  ,2.261576  , 

♦ 2.152733  ,2.038785  ,1.920001  ,1.796663  ,1.669369  ,1.537535  , 

♦ 1.902300  ,1.263769,1.122290  ,.9780996  ,.8315392  ,.6830565  . 

♦ . 5329537  , .38 15 867, .2293  197  , 7.699870  E-2/ 

DATA  W/l. 783276  E-3, 9. 196755  £-3,6.501866  E-3, 

♦8.839690  E-3, 1.112895  E-2, 1.336217  £-2,1.550370  E-2, 1.751906  E-2, 
♦1.939953  E-2, 2. 393731  E-2, 2. 222699  E- 2 , 2 . 3 1 398 3E -2 , 2 . 3601 36  E-2, 
♦2.353989  E-2, 2. 286831  E-2, 2. 153855  E- 2 , 1 . 99S706E- 2 , 1 .671 528  E-2, 
♦1.318996  E-2, 8. 997691  E-3, 9.098239  E- 3 , - 1 . 9 1 3 725  E-3 ,-7. 308989E-3, 

♦ -1. 397699  E-2, -1 .972312  E-2, - 2 . 585660  E- 2 ,- 3 . 1 69930  E -2  , -3 . 6879 39 
♦E-2,-9. 139236E-2,-9.936756  £-2,-9.730388  E-2 , -9 . 859856  E-2/ 

Y = 0, 

00  5 1=1,32 
Z=XII  I 

5 Y = y»ITlP»Z,TRANSI  ♦ T I P-Z,  TR  AN  SI  I ♦!«  1 II 
RETURN 
END 
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suiJPOuiiNE  080  coKPures  ^ ac-poitn  caussian  quaopature  over  the 

INTERVAL  AND  RETURNS  THE  ANSWER  Y.  SEE  06  EUR  THE  HEAMNG  OF  X ANO  W 
EXTERNAL  TRANS 
COKHCN/G/PI,PI2tZt<* 

0IKEN£!QN  XI^OItWIAOl 

DATA  - - ; X/3. 140191  ,3.134209  ,3.123458  ,3.10/950  , 

1 3.087710  ,3.062767  ,3.033161  ,2.998935  ,2.960143  ,2.916843  , 

2 2.869101  ,2.816939  ,2.760588  ,2.699983  ,2.635267  ,2.566537  , 

3 2.493899  ,2.417463  ,2.337346  ,2.253669  ,2.166561  ,2.0/6153  , 

4 1.982583  ,1.885994  ,1.786533  ,1.634352  ,1.5/9606  ,1.472454  , 

5 1.363063  ,1.251590  ,1.138214  ,1.323105  ,.9064376  ,.7383931  , 

6 .6691420  ,.5488749  ,.4277719  ,.3060175  ,.1837971  ,6.129682  E-2/ 

DATA  W/1  . 144949  E- 3, 2. 663461E-3  ,4 . 179626  E-3, 5. 687/02 

lE-3, 7. 132466  £-3,  8.656981  E-3,  1.010209  E-2, 1.150603  E-2, 1.285421  E 
2-2,1.412399  E-2, 1.530967  E-2, 1.637252  E-2, 1.729112  E-2 , 1 . 803665E-2 
3,1.857859  E-2,  1.888542  E-2, 1.892564  E-2, 1.866388  E-2, 1.808719  E-2, 
41.715643  E-2, 1.585776  E-2, 1.417906  E-2, 1.211641  E-2, 9. 675333  E-3, 6 
5.871675  E-3,  3.732713  E-3, 2. 971648  E-4 , -3 . 384 71 0 E-3 ,-7.25 1709  E-3, 
6-1.123276  E-2, -1.524848  E- 2, - 1 .92 1 330  E- 2 ,- 2 . 303786  E-2, 

7 -2.663188  E-2 ,- 2 . 990 720  E-2 ,- 3 . 2 78083  E-2 , -3 . 5 17804  E-2, 

8 -3.703520  E-2, -3 .830221  £-2,-3.894454  E-2/ 

Y = 3. 

DO  5 1=1,40 
Z = XI  I) 

5  y=Y*mPfZ,TRANS»  ♦ TlP-Z,TRANSn*wm 
RETURN 
END 


SUBROUTINE  C96I  Y,TRANSI 

SUBROUTINE  096  CCNPUTES  A- 96-POINT  GAUSSIAN  QUADRATURE  OVER-  THE 
INTERVAL  ANO  RETURNS  THE  ANSWER  Y.  SEE  06  FOR  THE  WEANING  OF  X ANO  W 
EXTERNAL  TRANS 
C0KM0N/C/PI,PI2,Z,P 
OIWENSION  X(46l,Wl4Si 

DATA  X/3. 143617  , 3. 1 3( 454  , 3 . 1 2 8969  ,3.118169  , 

1 3.104064  ,3.086669  ,3.066003  ,3.042088  ,3.314953  ,2.984616  , 

2 2.951119  ,2.914495  ,2.874782  ,2.832022  ,2.786261  ,2.737548  , 

3 2.685933  ,2.631472  ,2.574222  ,2.514244  ,2.451602  ,2.386361  , 

4 2.318592  ,2.248365  ,2.175756  ,2.1CC841  ,2.023699  ,1.944413  , 

5 1.863066  ,1.779745  ,1.694538  ,1.607535  ,1.518828  ,1.428512  , 

6 1.336682  ,1.243435  ,1.148871  ,1.C53C89  ,.9561907  ,.8582794  , 

7 .7594585  ,.6598328  ,.5595079  ,.4585899  ,.3571860  ,.2554036  , 

8 .1533505  ,5.113490  E-2/ 

DATA  M/7.967917  E-4, 1.853936  E-3, 2. 910500  E-3, 3. 963467 

• E-3,5. 010672  £-3,6.049410  E-3,7.076207  E-3, 8. 086678  E-3, 9. 075404 
♦E-3, 1.003584  E-2, 1.096023  E-2,1.18396C  E-2, 1.266371  E-2, 1.342112  E 
♦-2,1.409927  E-2, 1.468459  E-2, 1 .516267.  E-2, 1 . 55 1853  E-2, 1.573687  E- 
♦2,1.580239  E-2,t.57CC22  E-2, 1.541629  £-2,1.493764  E-2, 1.425388  E-2 
♦,1.335568  E-2, 1.223728  E-2, 1.089595  E-2, 9. 332676  E-3, 7. 552470  E-3, 
♦5.564738  E-3, 3.383468  E-3, 1.027341  E-3, -I .480296  E-3, -4. 111580  E-3 
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*,-6.03''i'i39  E-3, -9. 613338  E- 3 , - I . 2<i  0836  E - 2 , - 1 . 5 I ?89«  E-2  ,-l  . 788  169 

♦ C-2,-2.0'«728l  E-2,.-2.29CG8/  E -2 , - 2 . 5 1 <i  77  3 E-2  ,-2 . 7 7 3 E-2, -2. 88 

♦ 7858  £-2,-3  .030292  E- 2 , - 3 . I 39962  E- 2, - 3.2  I 3959  E-2  , -3 . 250807  E-2/ 

/=0. 

DO  5 1=1,98 
Z = XII  I 

5 Y=VU  I (P»7,TRANSI  ♦ T ( P-Z , TRANS  ))♦«(  II 
RETURN 
END 


